01 Neurodegenerative gene expression in peripheral immune cells

UCL

# Choose working location as jane or nemo
location <- "jane" 

if (location == "jane") {
  base_dir <- "/Volumes/lab-gandhis/home/users/wagena/periph_immune_neurodegen_publication/"
} else if (location == "nemo") {
  base_dir <- here::here()
  options(bitmapType='cairo') # Add graphics device for nemo
}

args <- list(PD_genes_file = str_c(base_dir, "/derived_data/Blauendraat_lancetneuro2020_table1_PDgenes.csv"),
             AD_genes_file = str_c(base_dir, "/derived_data/Goate_neurobiodis2020_AD_genes_w_destrooper.csv"),
             FTD_genes_file = str_c(base_dir, "/derived_data/Koch_ijms2023_FTD_genes_table1.csv"),
             aquino_eqtl_dir_path = file.path(base_dir, 
                                          "raw_data",
                                          "GWAS",
                                          "aquino_2023_covid_eqtl"),
             aquino_eqtl_summary_path =  file.path(base_dir, 
                                               "raw_data",
                                               "GWAS",
                                               "aquino_2023_covid_eqtl/aquino_summary_table.csv"),
             chen_2020_locus_to_gene_path = file.path(base_dir,
                                                      "raw_data",
                                                 "GWAS",
                                                 "chen_cell_2020",
                                                 "open_targets_locus_to_gene_neurodegen.csv"),
             chen_2020_substudy_lookup_table_path = file.path(base_dir,
                                                      "raw_data",
                                                 "GWAS",
                                                 "chen_cell_2020",
                                                 "chen_cell_2020_study_subset_lookup_table.tsv"),
             gwas_summary_numbers =  file.path(base_dir, "raw_data/gwas_summary_numbers_supp_figure.csv"),
             qtl_summary_numbers = file.path(base_dir, "raw_data/metadata_of_GWAS_QTL_included.csv"),
             gtex_gene_expression_per_tissue_file = file.path(base_dir, "raw_data/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct")
             
)



# Set defaults for ggplots 
theme_aw <- theme_bw(base_family = "Helvetica",
           base_size = 14) + 
    theme(panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          strip.background = element_blank(),
          legend.position = "bottom",
          axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          axis.title.y = element_text(vjust = 0.6),
          panel.spacing = unit(0.1, "lines"))


theme_set(theme_aw)

plot_labels = c("PD2019_ex23andMe" = "European PD GWAS",
                "foo2020" = "East Asian PD GWAS",
                 "rizig2023" = "African PD GWAS",
                "SNCA" = "SNCA",
                "GBA" = "GBA",
                "LRRK2" = "LRRK2",
                "BIN1" = "BIN1",
                "Monocyte_COV" = "Monocyte + SARS-CoV-2",
                "Monocyte_NS" = "Monocyte",
                "Dendritic-plasmacytoid_NS" = "Dendritic Plasmacytoid",
                "Mono-cd14_NS" = "Monocyte CD14+",
                "NS" = "-",
                "COV" = "SARS-CoV-2",
                "IAV" = "Influenza A",
                "T.CD4" = "CD4 T",
                "T.CD8" = "CD8 T",
                "NK" = "NK",
                "MYELOID" = "MYELOID",
                "B" = "B",
                "true_result" = "True result", 
                'all_gtex_genes' = "All genes", 
                "brain_gtex_genes" = "Brain genes")
                


plot_labeller = as_labeller(plot_labels)

    
my_palette <- c("Colocalised" = "darkslategrey",
             "Suggestive" = "aquamarine3",
             "Distinct" = "indianred3",
             "Not_signif" = "grey70",
             "PPH4" = "darkslategrey",
             "PPH3" = "indianred3",
             "Afr_As_Eur" = "indianred3",
             "As" = "darkslategrey",
             "Eur" = "aquamarine3",
             "Afr" = "deepskyblue3",
             "Afr_Eur" = "darkslateblue",
             "Expressed" = "darkslategrey",
             "Not Expressed" = "aquamarine3",
             "WBC association" = "darkslategrey",
             "Blood count association" = "darkslategrey",
             "No association" = "aquamarine3",
             "White blood cell count - negative β" = "slategrey",
              "White blood cell count - positive β" = "darkslategrey",
             "Monocyte count - negative β" =  "aquamarine3",
             "Monocyte count - positive β" =  "aquamarine4",
             "Lymphocyte counts - negative β" = "indianred3",
             "Lymphocyte counts - positive β" = "indianred4",
             "Neutrophil count - negative β" = "slateblue2",
              "Neutrophil count - positive β" = "darkslateblue",
             "Basophil count - negative β" = "slategrey",
             "Basophil count - positive β" = "darkslategrey",
             "Eosinophil counts - negative β" =  "aquamarine3",
             "Eosinophil counts - positive β" =  "aquamarine4",
             "Platelet count - negative β" = "indianred3",
             "Platelet count - positive β" = "indianred4",
             "Mean platelet volume - negative β" = "slateblue2",
              "Mean platelet volume - positive β" = "darkslateblue",
             "true_result" = "darkslategrey",
             'all_gtex_genes' = "indianred4", 
             "brain_gtex_genes" = "indianred3")

Aim: This analysis explores the expression of genes causally implicated in neurodegenerative diseases in peripheral immune cells.

1 Define genes of interest

1.1 Parkinson’s disease

The PD genes of interest come from Blaundraat Lancet Neurology 2020 table 1, using only the 14 genes that have high or very high confidence of being an actual PD gene (as defined by the authors).

PD_gene_df <- read_delim(args$PD_genes_file) %>% 
    dplyr::filter(Confidence_as_actual_PD_gene %in% c("Very high", "High")) %>% 
    dplyr::rename(Gene = `Empty Cell`) #%>% 
    # dplyr::mutate(Gene = dplyr::case_when(
    # Gene == "GBA" ~ "GBA1",  # Replace "GBA" with the proper name, "GBA1"
    # ))

 PD_genes <- tibble(disease = "PD",
                               gene = c(PD_gene_df$Gene, "RAB32")) # Remove KAT8 and CTSF as can't cherrypick here...
gt(PD_gene_df) %>% 
  opt_interactive()

1.2 Alzheimer’s disease

The AD genes of interest were derived from Alison Goates Genetic architecture of AD in Neurobiology of disease (2020). The genes selected were those highlighted in red, indicated increased confidence of causation.

AD_gene_df <- read_delim(args$AD_genes_file) %>% 
  dplyr::filter(Confident_goate == "Y") %>% 
  dplyr::select(-Confident_destrooper, -Reference)
AD_genes <- tibble(disease = "AD",
                   gene = AD_gene_df$`Reported gene`)

gt(AD_gene_df) %>% 
  opt_interactive(use_text_wrapping = T)

1.3 FTD genes

The FTD genes were derived from Antonioni et al Frontotemporal Dementia, Where Do We Stand? A Narrative Review. It takes the causative genes listed in table 1, and filters for those that account for at least 1% of disease.

FTD_gene_df <- read_delim(args$FTD_genes_file) %>% 
  dplyr::filter(Frequency %in% c("10%", "5%", "1%"))

FTD_genes <- tibble(disease = "FTD",
                   gene = FTD_gene_df$Gene) 

FTD_gene_df %>% 
  gt() %>% 
  opt_interactive(use_text_wrapping = T)

1.4 Total number of genes per disease

summary_genes <- bind_rows(PD_genes, AD_genes, FTD_genes) %>% 
  dplyr::mutate(disease = factor(disease, levels = c("AD", "PD", "FTD"))) %>% 
  arrange(disease, gene)
summary_gene_list <- summary_genes$gene %>% unlist()

write_csv(summary_genes,
          file = file.path(base_dir, "figures/summary_genes_table.csv"))

# # For thesis
# # Define each vector of genes
# AD_genes_vector <- AD_genes[["gene"]] %>% sort()
# PD_genes_vector <- PD_genes[["gene"]] %>% sort()
# FTD_genes_vector <- FTD_genes[["gene"]] %>% sort()
# 
# # Find the maximum length among the vectors
# max_length <- max(length(AD_genes_vector), length(PD_genes_vector), length(FTD_genes_vector))
# 
# # Pad each vector with NA to match the maximum length
# AD_genes_padded <- c(AD_genes_vector, rep(NA, max_length - length(AD_genes_vector))) 
# PD_genes_padded <- c(PD_genes_vector, rep(NA, max_length - length(PD_genes_vector))) 
# FTD_genes_padded <- c(FTD_genes_vector, rep(NA, max_length - length(FTD_genes_vector)))
# 
# # Combine the padded vectors into a tibble
# summary_genes_for_thesis <- tibble(
#   AD = AD_genes_padded,
#   PD = PD_genes_padded,
#   FTD = FTD_genes_padded
# ) 
# 
# summary_genes_for_thesis %>% 
#     write_csv(file = file.path(base_dir, "figures/summary_genes_for_thesis.csv"))
# 
# summary_genes_gt <- gt(summary_genes_for_thesis) %>% 
#     sub_missing(missing_text = "") %>% 
#     tab_style(
#     style = cell_text(weight = "bold"),
#     locations = cells_column_labels()
#     )
#     
# summary_genes_gt %>% 
#     gtsave(filename = file.path(base_dir, "figures/summary_genes_gt.png"))



summary_genes <- read_csv(file = file.path(base_dir, "figures/summary_genes_table.csv")) %>% 
  as_tibble()

summary_genes %>% 
  count(disease) %>% 
  dplyr::rename(`Number of genes` = n) %>% 
  gt()

disease	Number of genes
AD	18
FTD	6
PD	15

2 Studies included in analysis

2.1 eQTL Studies

4 different eQTL studies were included in the study. A multi-ancestry study (Aquino et al), an East Asian study (Ota et al), and 2 African/African admixed studies each covering a single peripheral immune cell: Nedelec et al exploring macrophages, and Quach et al exploring monocytes. All of these studies included peripheral immune cells that were stimulated, either in vitro using microbial or chemical stimuli, or, as in the case of Ota et al, in utilising human donors with a range of immune related diseases.

study_factor_order <- c("Aquino_2023",   "Ota_2021", "Nedelec_2016" ,  "Quach_2016")

study_summary <- read_csv(args$qtl_summary_numbers)

eqtl_summary <-  study_summary %>% 
  dplyr::filter(Study %in% c("Aquino_2023",   "Ota_2021", "Nedelec_2016" ,  "Quach_2016")) %>% 
  dplyr::mutate(Study = fct_relevel(Study, study_factor_order))


# Transform Stimulation into a factor with levels for better control in plotting
eqtl_summary$Ethnicities <- factor(eqtl_summary$Ethnicities)





# Separate the metrics into individual plots with specific aesthetics
individuals_plot <- ggplot(eqtl_summary, aes(x = "N_individuals", y = fct_rev(Study))) +
  geom_point(aes(size = N_individuals), colour = "darkslategrey") +  # Circle shape
  scale_size_continuous(range = c(1, 12), name = "N Individuals") +
  theme_minimal() +
  theme(axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.key.size =  unit(0.4, 'lines')) +
  labs(x = NULL, y = NULL)

celltypes_plot <- ggplot(eqtl_summary, aes(x = "N_celltypes", y = fct_rev(Study))) +
  geom_point(aes(size = N_celltypes), colour = "aquamarine3") +  # Square shape
  scale_size_continuous(range = c(1, 12), name = "N Celltypes") +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
         legend.key.size =  unit(0.4, 'lines')) +
  labs(x = NULL)


ethnicity_plot <- ggplot(eqtl_summary, aes(x = "Ancestry", y = fct_rev(Study))) +
  geom_point(aes(color = Ethnicities), shape = 15, size = 10) +  # Same shape for all ethnicities
  scale_color_manual(values = my_palette) +
  labs(x = NULL, y = NULL, color = "Ethnicities") +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          legend.key.size =  unit(0.4, 'lines'))


stimulation_plot <- ggplot(eqtl_summary, aes(x = "Stimulation", y = fct_rev(Study))) +
  geom_point(aes(size = ifelse(is.na(Stimulation), NA, 10)),  # Use conditional size
             shape = 15) +
  scale_size_continuous(range = c(0, 15), guide = "none") +  # Set scale for size, remove legend
  labs(x = NULL, y = NULL) +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          legend.key.size =  unit(0.4, 'lines'))

# Combine the plots
eqtl_combined_studies_plot <- (individuals_plot | celltypes_plot | ethnicity_plot | stimulation_plot) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom",
        legend.direction = "vertical")

# Print the combined plot
eqtl_combined_studies_plot

ggsave(eqtl_combined_studies_plot,
       file = file.path(base_dir, "figures/eqtl_datasets_summary.pdf"),
              units = "cm",
       height = 16,
       width = 12)

2.2 GWAS studies

gwas_studies <-  read_csv(args$gwas_summary_numbers) %>% 
  dplyr::mutate(Study = fct_relevel(Study,  c("Foo_2020",  "Nalls_2020", "Rizig_2023", "Jansen_2019", "Shigemizu_2021")))


# Transform Stimulation into a factor with levels for better control in plotting
gwas_studies$Ancestry <- factor(gwas_studies$Ancestry)



# Separate the metrics into individual plots with specific aesthetics
individuals_plot <- ggplot(gwas_studies, aes(x = "N_cases", y = fct_rev(Study))) +
  geom_point(aes(size = N_cases), colour = "darkslategrey") +  # Circle shape
  scale_size_continuous(range = c(0.5, 4), name = "N cases") +
  theme_minimal() +
  theme(axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.key.size =  unit(0.4, 'lines')) +
  labs(x = NULL, y = NULL)

celltypes_plot <- ggplot(gwas_studies, aes(x = "N_controls", y = fct_rev(Study))) +
  geom_point(aes(size = N_controls), colour = "aquamarine3") +  # Square shape
  scale_size_continuous(range = c(1, 16), name = "N controls") +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
         legend.key.size =  unit(0.4, 'lines')) +
  labs(x = NULL)


ethnicity_plot <- ggplot(gwas_studies, aes(x = "Ancestry", y = fct_rev(Study))) +
  geom_point(aes(color = Ancestry), shape = 15, size = 10) +  # Same shape for all ethnicities
  scale_color_manual(values = my_palette) +
  labs(x = NULL, y = NULL, color = "Ethnicities") +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          legend.key.size =  unit(0.4, 'lines'))


# Combine the plots
gwas_combined_studies_plot <- (individuals_plot | celltypes_plot | ethnicity_plot) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom",
        legend.direction = "vertical")

# Print the combined plot
gwas_combined_studies_plot

ggsave(gwas_combined_studies_plot,
       file = file.path(base_dir, "figures/gwas_datasets_summary.pdf"),
       units = "cm",
       height = 16,
       width = 12)

3 Which genes are expressed in peripheral immune cells: Aquino et al 2023

The Aquino et al study was chosen as the baseline study for analysis given that it incorporated donors from multiple ancestries, multiple cell types and activation states.

This is a summary figure from the original paper, showing the study structure and cell types explored.

Data utilised in this analysis was derived from:

• Supp Table 3c and 3d: gene expression counts (log2 count per million), of cell types (3c) and cell lineage (3d) after 6 hours simulation by virus
• Supp Table 4b showing variation in gene expression between African and European by cell lineage and activation states.

3.1 Table of genes expressed in peripheral immune cells

table_3d <- read_delim(str_c(args$aquino_eqtl_dir_path, "/aquino_supp_table_S3d.csv"),
                       skip = 3,
                       col_names = T)

# Update column names by combining first and second row
new_col_names <- paste(table_3d[1, ], table_3d[2, ], sep = "_")
names(table_3d) <- c("GeneID", "Symbol", new_col_names[-c(1,2)])  # Skip first two as they are GeneID and Symbol

# Remove the first two rows as they are now redundant
table_3d <- table_3d[-c(1, 2), ]

aquino_cell_order <- c("MYELOID",
                       "cDC",
                       "pDC",
                       "MONO CD14",
                       "MONO CD16",
                       "NK",
                       "NK CD56brt",
                       "NK CD56dim",
                       "NK MEM",
                       "ILC",
                       "T CD4",
                       "T CD4 E",
                       "T CD4 N",
                       "T CD8",
                       "T CD8 CM/EM",
                       "T CD8 EMRA",
                       "T CD8 N",
                       "Reg T",
                       "γδ T",
                       "MAIT",
                       "B",
                       "BMK",
                       "BML",
                       "BNK",
                       "BNL",
                       "Plasmablast")

# Convert the data to a long format
long_data <- table_3d %>%
  pivot_longer(
    cols = -c(GeneID, Symbol),  # Excluding GeneID and Symbol columns from the pivot
    names_to = c("Lineage", "virus"),
    names_sep = "_",  # Splitting based on the separator used in new column names
    values_to = "log2TPM"
  ) %>% 
  dplyr::filter(Symbol %in% summary_gene_list) %>% 
  dplyr::left_join(.,
                   summary_genes,
                   by = c("Symbol"="gene")) %>% 
  dplyr::mutate(log2TPM = as.numeric(log2TPM),
                virus = as_factor(virus) %>% fct_relevel(., c("NS", "IAV", "COV")),
                Lineage = as_factor(Lineage) %>%  fct_relevel(., aquino_cell_order)) %>% 
  group_by(Symbol) 

genes_in_aquino <- tibble(gene = long_data$Symbol) %>% 
  dplyr::distinct(gene) %>% 
  dplyr::mutate(gene_in_aquino = "Expressed")

                        
gene_expressed_in_data <- summary_genes %>% 
  dplyr::full_join(.,
            genes_in_aquino,
            by = c("gene")) %>% 
  dplyr::mutate(genes_in_aquino = case_when(gene_in_aquino == "Expressed" ~ "Expressed",
                                            TRUE ~ "Not Expressed")) %>% 
  dplyr::select(-gene_in_aquino)


gene_expressed_in_data %>% 
  gt() %>% 
  opt_interactive()

3.2 The instances of genes expressed per disease

aquino_gene_count_data <- gene_expressed_in_data %>% 
  count(disease, genes_in_aquino) 

aquino_gene_count <- aquino_gene_count_data %>% 
  ggplot(aes(x = n, y = fct_rev(disease), fill = fct_rev(genes_in_aquino))) +
  geom_bar(stat = "identity", position = "stack") +
  scale_fill_manual(values = my_palette,
                       guide = guide_legend(reverse = TRUE)) +
  labs(x = "Number of genes", y = NULL, fill = NULL) +
  theme_aw +
  theme(legend.position = "bottom", 
        panel.grid.major.x = element_line(),
        panel.grid.major.y = element_blank()
  )
        

aquino_gene_count

4 Association between peripheral blood counts and genetic variability at NDD genes: Chen et al 2020

To answer the question whether variants in NDD genes can influence Using Chen et al multi-ancestry phenotype QTL to explore whether peripheral blood counts are related to variation in NDD genes.

To do this, I have downloaded all associations from the locus-to-gene pipeline relating the the neurodegeneration genes above and this Chen paper via opentargets. Note, the metaanalysis incorporating the various ancestries do not have betas associated, just p values.

Out of the 32 genes that have expression in peripheral immune cells (based on aquino), I will look at which of those genes have associations with blood cell count.

4.1 Table of associations between blood count phenotype and NDD gene variants

traits_to_explore <- c("White blood cell count",
                       "Monocyte count",                 
                       "Lymphocyte counts",
                       "Neutrophil count",
                       "Basophil count",
                       "Eosinophil counts",
                       "Platelet count",
                       "Mean platelet volume"
                       )

opentargets <- read_delim(args$chen_2020_locus_to_gene_path) %>% 
    dplyr::left_join(.,
            gene_expressed_in_data %>% 
              dplyr::select(disease, gene,  genes_in_aquino),
            by = "gene") %>% 
  dplyr::filter(!gene %in% c("CTSB", "KAT8")) %>% 
  dplyr::select(study_id = `Study ID`,
                gene,
                disease,
                Trait = `Reported Trait`,
                N_participants = `Study N Initial`,  
                snp_ID = `Index Variant ID`,  
                snp_rs = `Index Variant RSID`,
                p_value = `P-Value`,
                beta = Beta,
                beta_ci_low = `Beta CI Lower`,
                beta_ci_high = `Beta CI Upper`,
                l2g_score = `L2G Pipeline Score`,
                has_sumstats = `Has sumstats`) %>% 
  dplyr::filter(Trait %in% traits_to_explore) %>% 
  dplyr::mutate(N_participants  = as.numeric(N_participants))



substudy_info <- read_delim(args$chen_2020_substudy_lookup_table_path) %>% 
  dplyr::select(study_id = accessionId,
                Trait = reportedTrait,
                discoverySampleAncestry,
                associationCount,
                summaryStatistics) %>% 
  extract(discoverySampleAncestry, into = c("N_participants", "Ancestry"), regex = "^([0-9]+)\\s(.+)$", remove = TRUE) %>% 
  dplyr::mutate(N_participants = as.numeric(N_participants)) %>% 
  arrange(study_id) %>% 
  dplyr::select(-Trait, -N_participants)

included_chen_study_ids <- substudy_info %>% 
    dplyr::filter(Ancestry %in% c("African American or Afro-Caribbean, African unspecified",
                                  "European",
                                  "East Asian",
                                  "Hispanic or Latin American",
                                  "South Asian",
                                  "African American or Afro-Caribbean, African unspecified, European, East Asian, Hispanic or Latin American, South Asian")) %>% 
    pull(study_id)

saveRDS(included_chen_study_ids,
        file = file.path(base_dir, "/derived_data/chen_derived_data/included_chen_study_ids.RDS"))


substudy_info$Ancestry %>%  unique()

## [1] "African American or Afro-Caribbean, African unspecified"                                                               
## [2] "European"                                                                                                              
## [3] "East Asian"                                                                                                            
## [4] "Hispanic or Latin American"                                                                                            
## [5] "South Asian"                                                                                                           
## [6] "African American or Afro-Caribbean, African unspecified, European, East Asian, Hispanic or Latin American, South Asian"

opentargets <- left_join(opentargets,
                           substudy_info,
                           by = c("study_id")) %>% 
  dplyr::relocate(study_id, Ancestry, N_participants, gene, disease, Trait, p_value, beta, beta_ci_low, beta_ci_high, everything()) %>% 
  dplyr::mutate(Ancestry = case_when(Ancestry == "African American or Afro-Caribbean, African unspecified, European, East Asian, Hispanic or Latin American, South Asian" ~ "Mixed",
                                     Ancestry == "East Asian" ~ "Asian",
                                     Ancestry == "European" ~ "European",
                                     Ancestry == "South Asian" ~ "South Asian",
                                     TRUE ~ "other")) %>% 
  dplyr::select(-summaryStatistics) %>% 
  arrange(Ancestry)

opentargets %>% 
 gt() %>% 
  opt_interactive(use_resizers = T,
                  use_filters = T
                  )

4.2 Plot of associations between NDD gene variants and peripheral white cell counts

mixed <- opentargets %>% dplyr::filter(Ancestry == "Mixed")
European <- opentargets %>%  dplyr::filter(Ancestry == "European")

# Create combined category of disease and direction of beta
opentargets <- opentargets %>% 
  dplyr::mutate(trait_direction = case_when(beta >0 ~ str_c(Trait, " - positive β"),
                                              beta<0 ~ str_c(Trait, " - negative β"),
                                              TRUE ~ NA_character_))

# # Check number of NDD genes that have an association in Chen in white cell count metrics
# opentargets %>% 
#   dplyr::filter(Trait %in% c("Monocyte count",                 
#                        "Lymphocyte counts",
#                        "Neutrophil count") ) %>% #,
#               #  Ancestry == "European") %>% 
#   dplyr::select(gene) %>% unique()


chen_association_count <- opentargets %>% 
  dplyr::filter(Trait %in% c("White blood cell count",
                       "Monocyte count",                 
                       "Lymphocyte counts",
                       "Neutrophil count")) %>% 
  ggplot(aes(x = gene, fill = trait_direction)) +
  geom_bar(stat = "count") +
  facet_nested( Trait + Ancestry ~ disease,
               scales = "free",
               space = "free") +
  scale_fill_manual(values = my_palette)+ 
  scale_y_continuous(breaks = function(x) seq(ceiling(x[1]), floor(x[2]), by = 1)) +
  theme_aw +
  theme(legend.position = "right",
        panel.grid.minor.y = element_blank()) +
  labs(fill = "Trait direction", y = "Count")


ggsave(chen_association_count,
       height = 34, width = 32, units = "cm",
       file = file.path(base_dir, "figures/chen_wbc_mixed_association_counts.pdf"))

ggsave(chen_association_count,
       height = 34, width = 32, units = "cm",
       file = file.path(base_dir, "figures/chen_wbc_mixed_association_counts.png"))


chen_association_count 

4.3 Instances where the NDD genes have an association with white cell count phenotypes.

# Create a tibble for genes by trait in Chen study
genes_in_chen <-  European %>% 
  count(gene, Trait) %>% 
  dplyr::mutate(chen_assoc = "Blood count association") 


# Defining the traits
traits <- c("White blood cell count", "Monocyte count", "Lymphocyte counts", "Neutrophil count")

# Create a data frame with all combinations of genes and traits
all_combinations <- expand.grid(gene = summary_genes$gene, Trait = traits)

# Join the all_combinations with the gene_expressed_in_chen table
expanded_table <- all_combinations %>% 
  dplyr::left_join(genes_in_chen, by = c("gene", "Trait")) %>% 
  dplyr::left_join(summary_genes, by = "gene") %>% 
  dplyr::mutate(genes_in_chen_relationship = ifelse(is.na(chen_assoc), "No association", chen_assoc))
            

chen_gene_count_data <- expanded_table %>% 
  group_by(disease) %>% 
  mutate(Trait = str_wrap(Trait, width = 16)) %>%
  count(disease, Trait, genes_in_chen_relationship)  %>% 
  dplyr::group_by(disease, Trait) %>%
  dplyr::summarise(
    total = sum(n),  # Sum of all entries for each group
    `Blood count association` = sum(n[genes_in_chen_relationship == "Blood count association"]),
    `No association` = sum(n[genes_in_chen_relationship == "No association"])
  ) %>%
  dplyr::mutate(
    percent_blood_count_association = `Blood count association` / total * 100,
    percent_no_association = `No association` / total * 100
  )


chen_gene_count_data %>%  gt()

Trait	total	Blood count association	No association	percent_blood_count_association	percent_no_association
AD
Lymphocyte counts	18	11	7	61.11111	38.88889
Monocyte count	18	10	8	55.55556	44.44444
Neutrophil count	18	12	6	66.66667	33.33333
White blood cell count	18	12	6	66.66667	33.33333
FTD
Lymphocyte counts	6	4	2	66.66667	33.33333
Monocyte count	6	3	3	50.00000	50.00000
Neutrophil count	6	3	3	50.00000	50.00000
White blood cell count	6	5	1	83.33333	16.66667
PD
Lymphocyte counts	15	4	11	26.66667	73.33333
Monocyte count	15	3	12	20.00000	80.00000
Neutrophil count	15	4	11	26.66667	73.33333
White blood cell count	15	5	10	33.33333	66.66667

chen_gene_count <- chen_gene_count_data %>% 
  dplyr::select(disease, Trait, `Blood count association`, `No association`) %>% 
  pivot_longer(cols = c("Blood count association", "No association"), names_to = "association", values_to = "count") %>% 
  ggplot(aes(x = disease, y = count, fill = fct_rev(association))) +
  geom_bar(stat = "identity", position = "fill") +
  scale_fill_manual(values = my_palette) +
  labs(x = NULL, y = "Proportion of genes", fill = NULL) +
  theme_aw +
  theme(legend.position = "bottom") +
  facet_grid(~ Trait)

chen_gene_count

5 How are NDD genes expressed across cell types and activation states: Aquino et al.

5.1 Proportion of gene expression for each gene, normalised by cell type & activation state

This section explores the pattern of expression of NDD genes in peripheral immune cells. To identify which peripheral cell has the greatest expression, the results were normalised to the cell type and activation state with the greatest expression per gene. The figure shows the activation state with the maximum proportion per cell type.

table_3c <- read_delim(str_c(args$aquino_eqtl_dir_path, "/aquino_supp_table_S3c.csv"),
                       skip = 3,
                       col_names = T)

# Update column names by combining first and second row
new_col_names <- paste(table_3c[1, ], table_3c[2, ], sep = "_")
names(table_3c) <- c("GeneID", "Symbol", new_col_names[-c(1,2)])  # Skip first two as they are GeneID and Symbol

# Remove the first two rows as they are now redundant
table_3c <- table_3c[-c(1, 2), ]

aquino_cell_order <- tibble(Cell = c("cDC",
                                     "pDC",
                                     "MONO CD14",
                                     "MONO CD16",
                                     "NK CD56brt",
                                     "NK CD56dim",
                                     "NK MEM",
                                     "ILC",
                                     "T CD4 E",
                                     "T CD4 N",
                                     "Reg T",
                                     "T CD8 CM/EM",
                                     "T CD8 EMRA",
                                     "T CD8 N",
                                     "MAIT",
                                     "γδ T",
                                     "BMK",
                                     "BML",
                                     "BNK",
                                     "BNL",
                                     "Plasmablast"),
                            Lineage = c("Dendritic",
                                        "Dendritic",
                                        "Monocyte",
                                        "Monocyte",
                                        "NK",
                                        "NK",
                                        "NK",
                                        "NK",
                                        "CD4",
                                        "CD4",
                                        "CD4",
                                        "CD8",
                                        "CD8",
                                        "CD8",
                                        "CD8",
                                        "CD8",
                                        "B",
                                        "B",
                                        "B",
                                        "B",
                                        "B")
)


# Convert the data to a long format
long_data <- table_3c %>%
    pivot_longer(
        cols = -c(GeneID, Symbol),  # Excluding GeneID and Symbol columns from the pivot
        names_to = c("Cell", "virus"),
        names_sep = "_",  # Splitting based on the separator used in new column names
        values_to = "log2TPM" 
    ) %>% 
    dplyr::mutate(log2TPM = as.numeric(log2TPM),
                  TPM = 2^log2TPM)  # Reversing the log2 transformation to get the original TPM values


filtered_data <- long_data %>%
  dplyr::filter(Symbol %in% summary_gene_list) %>%
  dplyr::left_join(.,
                   summary_genes,
                   by = c("Symbol"="gene")) %>%
  dplyr::left_join(.,
                   aquino_cell_order,
                   by = "Cell") %>%
  dplyr::mutate(#log2TPM = as.numeric(log2TPM),
                virus = fct_relevel(virus, c("NS", "IAV", "COV")),
                Cell = fct_relevel(Cell, aquino_cell_order$Cell),
                Lineage = fct_relevel(Lineage, aquino_cell_order$Lineage %>% unique())) %>%
  group_by(Symbol) %>%
  mutate(
   #  TPM = 2^log2TPM,  # Reversing the log2 transformation to get the original TPM values
    proportion_expression_per_gene = TPM / max(TPM)) %>%   # Calculating the proportion to allow comparison within a gene
  ungroup()

# Function to wrap text at 20 characters
wrap_text <- function(x) {
  sapply(x, function(y) str_wrap(y, width = 8))  # Adjust width as needed
}



# Legend width paremeter
legend_width <- unit(1.5, "lines")

max_proportion_per_treatment_data <- filtered_data %>% 
  as_tibble() %>% 
    dplyr::group_by(Symbol, Cell) %>% 
  dplyr::slice_max(proportion_expression_per_gene, n = 1, with_ties = FALSE) %>%   # Select max proportional value out of the NS, COV and IAV treatments for each celltype and gene
  ungroup() 

dendritic_cell_genes <- c("APP", "PSEN2", "PARK7")

# Find ratio of values between CD14 monocytes and CD56dim NK cells (for whichever treatment), to explore the pattern between different genes
mono_nk_ratio <- max_proportion_per_treatment_data %>%
  dplyr::filter(Cell %in% c("MONO CD14", "NK CD56dim")) %>% 
  dplyr::select(-Lineage, -TPM, -log2TPM, -virus) %>% 
  pivot_wider(names_from = Cell,
              values_from = proportion_expression_per_gene) %>% 
  dplyr::mutate(log_mono_nk_ratio = log10(`MONO CD14` / `NK CD56dim`)) %>% 
  arrange(log_mono_nk_ratio) %>% 
  dplyr::select(Symbol, disease, log_mono_nk_ratio, cd14_mono_prop = `MONO CD14`, cd56dim_nk_prop = `NK CD56dim`)
  

max_proportion_per_treatment_data <- max_proportion_per_treatment_data %>% 
  left_join(.,
            mono_nk_ratio,
            by = c("Symbol", "disease")) %>% 
  dplyr::mutate(Symbol = fct_reorder(Symbol, -cd14_mono_prop), # Reorder factors by descencing CD14 Monocyte proportion values
                Symbol = fct_relevel(Symbol, dendritic_cell_genes, after = 0))


cell_factor_levels <- levels(max_proportion_per_treatment_data$Symbol)


aquino_cell_expression_plot_max <- max_proportion_per_treatment_data  %>% 
  ggplot(aes(x = Cell, y = fct_rev(Symbol), fill = proportion_expression_per_gene)) +
  geom_tile() +
  facet_nested(disease ~ Lineage,
             scales = "free",
             space = "free",
             labeller = labeller(Lineage = wrap_text),
              ) +  # Applying wrap_text to Lineage labels
  theme_aw +
  theme(legend.key.size = legend_width,
        legend.position = "bottom",
          strip.background.x  = element_rect(color="black", linewidth=1, linetype="solid"),
        legend.title = element_text(vjust = 0.7)) +  # Center-align legend title
  scale_fill_viridis_c() +
  labs(title = "Proportion of gene expression per gene",
       subtitle = "Maximal proportion over treatments, ordered by CD14-Mono and DC", 
       fill = "Proportion expression",
       x = NULL,
       y = "Gene")


aquino_cell_expression_plot_max

5.2 Enrichment of expression of NDD genes in peripheral cells, relative to random/brain related genes

To explore whether NDD gene expression is enriched in peripheral immune cells (the Aquino dataset containing 12676 unique genes), relative to either a random set of all gtex genes (56200 genes), or cortical brain genes (29560 genes). This uses a summary of gene expression (median across samples) of tissues from GTEXv8.

gtex_gene_expression <- read_delim(args$gtex_gene_expression_per_tissue_file,
                                   skip =2) %>% 
    dplyr::select(gene = Name,
                  symbol = Description,
                  brain_cortex = `Brain - Cortex`) %>% 
    dplyr::mutate(gene = str_remove(gene, "\\.[0-9]+$"),  # Remove the suffix from the gene ID from gtex
                  symbol = str_trim(symbol),
                  symbol = dplyr::case_when(
                      symbol == "GBA1" ~ "GBA",
                      symbol == "PARK2" ~ "PRKN",
                      TRUE ~ symbol))  # Replace "GBA1" with the gene name used in Aquino, "GBA"


# 1) Set up target genes
# Get ENSG numbers of the 39 NDD genes from GTEX dataset
NDD_geneID <- gtex_gene_expression %>% 
    dplyr::filter(symbol %in% summary_gene_list) %>% 
    dplyr::pull(gene) %>% sort()

saveRDS(NDD_geneID,
        file = file.path(base_dir, "/derived_data/NDD_ENSG_gene_ids.RDS"))

NDD_geneID_in_aquino <- gtex_gene_expression %>% 
    dplyr::filter(symbol %in% genes_in_aquino$gene) %>% 
    dplyr::pull(gene) %>% sort()

saveRDS(NDD_geneID_in_aquino,
        file = file.path(base_dir, "/derived_data/NDD_ENSG_gene_ids_in_aquino.RDS"))
        
# Number of target genes
target_genes <- NDD_geneID
n_genes_tested <- length(target_genes)

# Identify all the expressed aquino genes and number expressed
expressed_target_genes <- genes_in_aquino$gene 
n_aquino_target_genes_expressed <- length(expressed_target_genes)



# 2) Set up aqiuno dataset
aquino_data_filtered <- long_data %>% 
     dplyr::select(GeneID, Symbol, Cell, virus, TPM)
expressed_aquino_genes <- aquino_data_filtered$GeneID %>%  unique()
m <- length(expressed_aquino_genes)

# Set up backgrounds from GTEX
# Create list of all 56200 gtex genes
all_gtex_genes <- gtex_gene_expression %>% 
    dplyr::pull(gene)

saveRDS(all_gtex_genes,
        file = file.path(base_dir, "/derived_data/gtex_all_ENSG_ids.RDS"))

# Create list of 29560 brain gtex genes
brain_gtex_genes <- gtex_gene_expression %>% 
    dplyr::filter(brain_cortex >0) %>% 
    dplyr::pull(gene)

saveRDS(brain_gtex_genes,
        file = file.path(base_dir, "/derived_data/gtex_brain_ENSG_ids.RDS"))

enrichment_backgrounds <- list(all_gtex_genes = all_gtex_genes,
                               brain_gtex_genes = brain_gtex_genes)




# Initialize list to store results for each background
results_list <- list()

n_iterations = 1000

 for (bg_name in names(enrichment_backgrounds)) {
   # bg_name <- names(enrichment_backgrounds)[1]
    gene_background <- enrichment_backgrounds[[bg_name]]
    
    # Background genes excluding target genes
    background_only_genes <- setdiff(gene_background, target_genes)
    n <- length(background_only_genes)
    
    # Set seed for reproducibility
    set.seed(613)
    
    # Perform random sampling using replicate for efficiency
    random_counts <- replicate(n_iterations, {
      # Randomly sample genes from the background
      random_sample <- sample(background_only_genes, size = n_genes_tested, replace = FALSE)
      
      # Count the number of sampled genes that are expressed
      sum(random_sample %in% expressed_aquino_genes)
    })
    
    # Median number of expressed genes in random samples
    median_random_expressed <- median(random_counts)
    
    # Standard deviation
    sd_random_expressed <- sd(random_counts)
    
    # 95% Confidence Interval
    ci_lower <- quantile(random_counts, 0.025)
    ci_upper <- quantile(random_counts, 0.975)
    
    # Empirical p-value
    p_value_empirical <- sum(random_counts >= n_aquino_target_genes_expressed) / n_iterations
    
    # Hypergeometric test p-value
    # Of the format: phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE), where:
    #   q = number of target genes (successes) in a bootstrap sample, 
    #   m is number of target genes in the background population,
    #   n is the total size of the background population
    #   k is the size of the sample
    
    k <- n_genes_tested # n is the size of the sample drawn
    q <- n_aquino_target_genes_expressed # The number of successes in the current bootstrap sample
    p_value_hyper <- phyper(q, m, n, k, lower.tail = FALSE)
    
    # Compile results into a dataframe
    result_df <- tibble(
        background = bg_name,
        n_genes_tested = n_genes_tested,
        n_target_genes_expressed = n_aquino_target_genes_expressed,
        median_random_expressed = median_random_expressed,
        sd_random_expressed = sd_random_expressed,
        ci_lower = as.numeric(ci_lower),
        ci_upper = as.numeric(ci_upper),
        p_value_empirical = p_value_empirical,
        p_value_hyper = p_value_hyper,
        raw_values = list(random_counts)
    )
    
    results_list[[bg_name]] <- result_df 
    
  }

results_df <- do.call(rbind, results_list)

gt(results_df)

background	n_genes_tested	n_target_genes_expressed	median_random_expressed	sd_random_expressed	ci_lower	ci_upper	p_value_empirical	p_value_hyper	raw_values
all_gtex_genes	39	32	8	2.608319	4	14	0	5.463316e-19	11, 7, 10, 11, 8, 8, 10, 6, 14, 9, 7, 10, 12, 2, 6, 8, 7, 8, 4, 10, 10, 4, 6, 5, 11, 7, 6, 8, 8, 8, 11, 8, 8, 6, 11, 4, 9, 6, 9, 8, 5, 9, 6, 12, 12, 5, 6, 7, 9, 7, 7, 8, 8, 6, 10, 8, 10, 10, 10, 9, 10, 13, 8, 7, 14, 13, 9, 5, 13, 9, 9, 14, 9, 12, 11, 11, 12, 7, 8, 4, 6, 8, 11, 7, 11, 12, 10, 9, 8, 7, 8, 4, 10, 10, 8, 11, 10, 5, 11, 9, 7, 6, 6, 7, 7, 8, 9, 8, 7, 8, 4, 9, 6, 13, 11, 14, 4, 5, 9, 7, 10, 12, 9, 12, 6, 15, 8, 7, 10, 11, 8, 9, 6, 10, 5, 9, 12, 15, 6, 10, 13, 11, 8, 12, 10, 11, 9, 10, 11, 9, 8, 7, 7, 9, 6, 5, 9, 11, 12, 8, 6, 6, 8, 7, 6, 8, 6, 10, 14, 7, 6, 11, 12, 11, 13, 9, 11, 8, 7, 7, 8, 10, 8, 8, 12, 7, 5, 9, 9, 8, 11, 9, 10, 7, 11, 8, 12, 9, 13, 8, 6, 7, 11, 5, 9, 6, 9, 12, 10, 13, 11, 7, 9, 11, 5, 6, 8, 8, 9, 7, 8, 9, 13, 8, 5, 8, 10, 9, 8, 7, 8, 7, 10, 9, 4, 3, 9, 10, 11, 9, 8, 10, 8, 12, 9, 12, 12, 8, 9, 5, 7, 12, 10, 11, 11, 9, 3, 12, 9, 9, 11, 7, 9, 8, 6, 11, 5, 7, 9, 7, 8, 9, 10, 9, 5, 8, 10, 14, 10, 9, 9, 10, 8, 15, 5, 6, 9, 12, 6, 8, 12, 9, 13, 7, 8, 9, 7, 10, 9, 9, 3, 4, 4, 4, 4, 10, 8, 10, 7, 4, 6, 9, 5, 8, 12, 8, 10, 7, 4, 9, 10, 7, 9, 5, 6, 12, 15, 11, 7, 3, 9, 7, 12, 11, 11, 11, 11, 9, 11, 10, 6, 7, 11, 7, 14, 6, 12, 14, 12, 6, 4, 8, 7, 8, 11, 5, 13, 11, 6, 8, 7, 9, 9, 11, 9, 9, 9, 8, 9, 3, 12, 13, 6, 11, 9, 4, 9, 9, 9, 8, 9, 16, 11, 12, 6, 9, 6, 10, 8, 6, 7, 10, 9, 10, 8, 11, 7, 10, 10, 10, 8, 10, 10, 8, 7, 9, 13, 8, 7, 7, 4, 9, 10, 11, 7, 9, 2, 7, 8, 6, 8, 7, 8, 6, 8, 8, 10, 7, 13, 8, 14, 11, 8, 11, 8, 10, 7, 11, 7, 8, 11, 13, 9, 11, 4, 11, 9, 7, 6, 9, 7, 8, 12, 11, 4, 8, 5, 9, 14, 5, 8, 8, 12, 8, 7, 13, 7, 7, 6, 9, 8, 17, 9, 12, 6, 10, 10, 10, 11, 5, 9, 9, 8, 11, 5, 9, 6, 13, 11, 9, 8, 8, 7, 7, 6, 12, 6, 10, 11, 11, 8, 6, 16, 6, 11, 11, 6, 10, 7, 4, 9, 7, 7, 14, 10, 15, 8, 6, 4, 7, 9, 9, 8, 9, 2, 12, 7, 5, 5, 6, 11, 11, 5, 12, 9, 2, 13, 9, 9, 10, 11, 6, 10, 11, 9, 6, 5, 11, 9, 9, 11, 12, 9, 9, 6, 8, 10, 5, 9, 9, 4, 10, 13, 8, 4, 8, 14, 9, 6, 5, 9, 6, 5, 8, 4, 11, 14, 7, 7, 8, 10, 6, 6, 6, 5, 15, 13, 9, 7, 9, 10, 5, 10, 6, 11, 7, 7, 7, 11, 17, 10, 9, 8, 6, 5, 10, 9, 12, 6, 8, 7, 5, 15, 8, 11, 10, 11, 8, 7, 8, 5, 8, 8, 9, 10, 5, 6, 13, 8, 11, 6, 10, 8, 8, 6, 3, 8, 9, 9, 6, 9, 9, 14, 12, 10, 9, 3, 7, 6, 8, 4, 14, 7, 7, 5, 11, 7, 17, 10, 9, 13, 15, 8, 8, 8, 7, 9, 8, 12, 10, 5, 8, 7, 7, 6, 7, 11, 10, 12, 9, 8, 6, 9, 7, 9, 6, 8, 11, 10, 6, 7, 4, 14, 5, 9, 7, 4, 3, 7, 11, 8, 7, 15, 5, 10, 7, 4, 10, 8, 10, 11, 2, 9, 13, 6, 8, 10, 5, 4, 6, 9, 6, 8, 7, 7, 9, 8, 13, 8, 5, 11, 12, 10, 5, 6, 6, 10, 11, 11, 9, 9, 13, 8, 6, 6, 6, 12, 10, 7, 7, 11, 8, 5, 10, 10, 11, 9, 15, 11, 7, 9, 10, 9, 5, 6, 6, 13, 12, 7, 9, 7, 8, 4, 6, 10, 12, 10, 8, 8, 8, 7, 7, 12, 13, 4, 9, 8, 12, 10, 9, 10, 7, 8, 11, 10, 9, 5, 5, 8, 9, 6, 12, 10, 6, 12, 7, 10, 9, 8, 8, 7, 10, 6, 15, 12, 12, 11, 8, 11, 5, 9, 9, 6, 12, 8, 14, 6, 6, 11, 7, 8, 10, 7, 8, 7, 8, 7, 13, 6, 7, 5, 7, 10, 6, 10, 9, 3, 7, 8, 6, 7, 13, 8, 5, 12, 8, 10, 6, 11, 8, 6, 8, 9, 12, 9, 7, 10, 9, 7, 13, 8, 8, 6, 13, 8, 4, 10, 4, 11, 6, 3, 7, 11, 10, 11, 5, 12, 8, 10, 13, 12, 11, 8, 7, 8, 6, 11, 12, 9, 7, 6, 7, 8, 10, 9, 8, 5, 11, 7, 8, 10, 8, 11, 9, 9, 3, 10, 8, 13, 11, 6, 5, 4, 10, 13, 8, 10, 12, 7, 9, 11, 8, 13, 4, 9, 8, 7, 12, 7, 7, 7, 10, 9, 6, 8, 9, 7, 6, 13, 4, 10, 4, 8, 6, 8, 13, 12, 9, 5, 6, 11, 6, 9, 4, 11, 5, 8, 6, 4, 7, 10, 9, 8, 7, 5, 9, 7, 10, 7, 4, 10, 7, 7, 8, 12, 10, 9, 9, 8, 5, 13, 12, 9, 12, 6, 6, 8, 11, 10, 13
brain_gtex_genes	39	32	16	3.021630	10	22	0	2.342556e-12	14, 17, 11, 11, 12, 15, 17, 15, 22, 19, 15, 16, 15, 13, 16, 14, 16, 19, 18, 15, 14, 12, 21, 14, 20, 21, 19, 14, 14, 19, 16, 13, 14, 16, 12, 18, 16, 17, 14, 20, 16, 14, 14, 14, 16, 16, 11, 16, 14, 17, 17, 18, 15, 23, 18, 16, 11, 19, 19, 13, 17, 18, 19, 17, 19, 16, 13, 11, 22, 12, 14, 13, 17, 21, 13, 21, 10, 15, 15, 11, 13, 14, 16, 18, 17, 11, 15, 12, 9, 19, 17, 18, 21, 16, 19, 13, 11, 16, 16, 18, 14, 13, 16, 11, 10, 17, 13, 19, 18, 12, 19, 12, 19, 16, 19, 11, 16, 15, 22, 20, 14, 13, 17, 22, 20, 12, 18, 15, 16, 15, 16, 16, 14, 15, 19, 14, 18, 15, 15, 15, 8, 18, 18, 16, 16, 14, 23, 13, 16, 11, 17, 20, 19, 9, 21, 24, 15, 15, 10, 13, 15, 15, 19, 15, 13, 17, 17, 16, 18, 12, 12, 13, 21, 19, 18, 18, 14, 26, 16, 14, 17, 17, 18, 15, 16, 15, 17, 14, 19, 19, 16, 15, 16, 17, 20, 17, 11, 16, 20, 13, 19, 13, 17, 23, 11, 15, 16, 17, 16, 12, 21, 17, 20, 21, 16, 15, 13, 11, 15, 21, 21, 17, 15, 19, 13, 12, 14, 17, 14, 16, 18, 14, 14, 13, 12, 23, 12, 21, 13, 12, 18, 22, 16, 18, 14, 11, 15, 15, 16, 15, 15, 19, 17, 15, 14, 14, 14, 19, 20, 17, 16, 15, 17, 16, 18, 16, 16, 15, 15, 22, 15, 11, 16, 19, 12, 13, 15, 15, 22, 20, 15, 14, 14, 14, 19, 17, 24, 17, 13, 16, 12, 18, 16, 12, 15, 15, 18, 16, 9, 16, 18, 18, 17, 15, 17, 16, 19, 13, 17, 23, 13, 19, 24, 9, 19, 11, 15, 15, 17, 11, 21, 17, 14, 20, 19, 15, 20, 15, 15, 19, 13, 14, 21, 16, 18, 18, 14, 13, 12, 17, 16, 17, 19, 9, 16, 19, 16, 14, 19, 9, 19, 18, 21, 15, 21, 22, 14, 13, 19, 16, 20, 10, 17, 18, 17, 17, 17, 19, 17, 13, 15, 16, 11, 16, 18, 15, 12, 13, 16, 19, 22, 19, 13, 23, 15, 19, 20, 19, 17, 17, 17, 17, 14, 9, 18, 18, 14, 20, 18, 16, 19, 20, 18, 21, 15, 14, 17, 17, 12, 13, 14, 13, 17, 20, 16, 17, 14, 18, 11, 16, 14, 19, 12, 12, 14, 16, 20, 15, 14, 17, 18, 13, 14, 14, 10, 13, 18, 18, 21, 23, 13, 18, 15, 14, 17, 13, 18, 14, 18, 19, 15, 17, 13, 16, 19, 13, 17, 15, 10, 18, 18, 18, 15, 16, 20, 17, 14, 14, 18, 10, 15, 16, 17, 13, 12, 25, 17, 16, 17, 20, 16, 15, 21, 14, 18, 19, 10, 19, 12, 15, 19, 18, 16, 16, 14, 15, 17, 10, 17, 16, 16, 15, 13, 17, 23, 16, 19, 16, 20, 14, 15, 17, 14, 13, 10, 17, 17, 11, 18, 19, 19, 12, 19, 22, 19, 19, 12, 19, 22, 18, 14, 18, 11, 18, 20, 14, 20, 16, 18, 15, 12, 19, 19, 18, 21, 19, 17, 17, 13, 19, 18, 16, 17, 15, 12, 18, 12, 11, 15, 16, 16, 16, 10, 15, 18, 15, 19, 14, 20, 15, 21, 14, 17, 17, 16, 10, 17, 17, 17, 15, 12, 13, 15, 16, 22, 20, 16, 20, 16, 15, 19, 13, 24, 18, 20, 12, 17, 17, 16, 17, 12, 18, 19, 14, 12, 19, 13, 16, 17, 18, 14, 14, 16, 18, 17, 14, 13, 17, 12, 16, 19, 16, 16, 14, 17, 19, 13, 13, 17, 18, 19, 18, 20, 15, 16, 17, 16, 15, 13, 16, 16, 18, 19, 19, 15, 14, 17, 15, 17, 14, 13, 13, 21, 9, 15, 20, 14, 15, 16, 20, 23, 14, 19, 12, 15, 23, 15, 22, 22, 9, 18, 15, 14, 19, 19, 18, 17, 20, 19, 17, 17, 13, 14, 19, 19, 15, 15, 14, 17, 17, 15, 18, 17, 20, 21, 16, 15, 14, 16, 11, 13, 20, 16, 12, 10, 14, 14, 15, 14, 14, 12, 15, 18, 17, 10, 24, 15, 18, 17, 15, 15, 22, 13, 18, 13, 13, 14, 18, 14, 15, 16, 15, 15, 17, 18, 13, 16, 15, 16, 20, 18, 21, 15, 16, 20, 14, 14, 16, 13, 21, 19, 18, 12, 17, 18, 14, 14, 16, 19, 19, 18, 15, 18, 15, 22, 14, 13, 13, 23, 15, 7, 18, 18, 17, 14, 17, 13, 22, 12, 18, 14, 13, 17, 18, 20, 20, 17, 18, 21, 15, 18, 15, 14, 18, 21, 12, 11, 15, 14, 17, 10, 11, 18, 17, 21, 12, 13, 13, 18, 17, 16, 18, 19, 13, 12, 24, 15, 17, 19, 20, 15, 14, 20, 16, 17, 14, 17, 18, 14, 19, 15, 12, 10, 24, 15, 11, 18, 17, 16, 12, 15, 15, 20, 17, 18, 15, 17, 19, 14, 20, 14, 21, 15, 22, 18, 19, 16, 13, 18, 12, 19, 17, 18, 18, 13, 14, 20, 16, 17, 15, 17, 13, 13, 15, 16, 17, 13, 13, 13, 18, 15, 16, 13, 18, 18, 20, 14, 18, 15, 12, 19, 13, 12, 16, 16, 14, 19, 19, 14, 22, 19, 20, 16, 18, 18, 19, 13, 15, 13, 12, 18, 18, 13, 12, 18, 15, 14, 15, 16, 16, 21, 14, 17, 16, 13, 14, 12, 17, 15, 19, 16, 16, 17, 18, 16, 17, 21, 20, 17, 22, 18, 17, 20, 15, 20, 20, 17, 13, 13, 15, 12, 18, 16, 20, 12, 17, 15, 21, 18, 16, 19, 15, 15, 13, 14, 20, 19, 13, 16, 15, 16, 16, 15, 18, 18, 12, 13, 16, 17, 12, 21, 10, 16, 16, 11, 15, 15, 8, 15, 20, 19, 17, 21, 18, 19, 14, 19, 15, 21, 21

# Summary dataframe for plotting
# true_results <- tibble(background = "true_result",
#               median_expression = n_aquino_target_genes_expressed,
#               ci_lower = n_aquino_target_genes_expressed,
#               ci_upper = n_aquino_target_genes_expressed)
# 
# enrichment_plot_data <- results_df %>% 
#     dplyr::select(background, median_expression = median_random_expressed, ci_lower, ci_upper) %>% 
#     bind_rows(true_results) %>% 
#     dplyr::mutate(background = fct_relevel(background, c("true_result", 'all_gtex_genes', "brain_gtex_genes"))) 

# enrichment_plot <- ggplot(enrichment_plot_data) +
#     geom_pointrange(aes(x = "", y = median_expression, ymin = ci_lower, ymax = ci_upper, colour = background), 
#                     position = position_dodge(width = 0.1)) +
#     scale_colour_manual(values = my_palette,  labels = plot_labeller) +
#     labs(x = NULL, y = "PIC-expressed genes", colour = NULL) +
#     scale_y_continuous(breaks = c(10,20,30,39), limits = c(0,37.2)) +
#     theme(axis.ticks.x = element_blank(),
#           legend.direction = "vertical")
# 
# enrichment_plot

enrichment_histogram_data <- results_df %>% 
    unnest(cols = c(raw_values)) %>% 
    dplyr::select(background, raw_values) 

median_data <- enrichment_histogram_data %>%
  group_by(background) %>%
  summarize(
    median_value = median(raw_values),
    max_density = max(density(raw_values)$y)
  )

   
enrichment_expression <- enrichment_histogram_data %>% 
    ggplot(aes(x = raw_values, colour = background)) + 
    geom_density(linewidth = 1.5, show.legend = FALSE) +
    geom_vline(xintercept = 32, colour = "darkslategrey", linewidth = 1, linetype = "dashed") +
    annotate("text", x = 32, y = Inf, size = 5, label = "\nExpressed genes", hjust = 1.1, angle= 90, colour = "darkslategrey") +
    annotate("text", x = 16, y = 0.15, 
             label = "p-value < 0.001", 
             colour = "indianred4", size = 5) +
    annotate("text", x = 24, y = 0.12,
             label = "p-value < 0.001", 
             colour = "indianred3", size = 5) +
    scale_x_continuous(limits = c(0,39),
                       breaks = c(0,10, 20, 30, 39)) +
    scale_colour_manual(values = my_palette, labels = plot_labeller) +
    labs(x = "Number of genes", y = "Density", colour = "Background") +
    geom_segment(data = median_data, aes(x = median_value, 
                                         xend =  median_value, 
                                         colour = background,
                                         y = 0, yend = max_density), 
                 linetype = "dashed", linewidth = 1) +
    guides(colour = guide_legend(override.aes = list(linetype = "solid", size = 1.5))) +
    theme(legend.key = element_rect(fill = NA, color = NA),
          panel.grid.major.x = element_line(),
          panel.grid.major.y = element_blank(),
          panel.grid.minor.y = element_blank()
    )





total_gene_counts <- gene_expressed_in_data %>% 
    count(genes_in_aquino) %>% 
    dplyr::mutate(disease = "Total")

total_aquino_gene_count <- total_gene_counts %>% 
    ggplot(aes(x = n, y = disease, fill = fct_rev(genes_in_aquino))) +
    geom_bar(stat = "identity", position = "stack") +
    scale_fill_manual(values = my_palette,
                      guide = guide_legend(reverse = TRUE)) +
    scale_x_continuous(limits = c(0,39),
                       breaks = c(0, 10, 20, 30, 39)) +
labs(x = "Number of genes", y = NULL, fill = NULL) +
    theme_aw +
    theme(legend.position = "bottom", 
          panel.grid.major.x = element_line(),
          panel.grid.major.y = element_blank()
    )

enrichment_plot <- (aquino_gene_count | 
                  (
                    (
                      enrichment_expression + 
                      labs(x = NULL) +
                      theme(axis.text.x = element_blank())
                    ) / 
                    total_aquino_gene_count
                  ) + 
                  plot_layout(heights = c(6, 1))
              ) + 
              plot_layout(guides = "collect") +
    plot_annotation(tag_levels = c("a", "b"))

enrichment_plot

ggsave(enrichment_plot,
       file = file.path(base_dir, "figures/thesis_gene_expression_enrichment.pdf"),
       height = 14, width = 30, units = "cm")

ggsave(enrichment_plot,
       file = file.path(base_dir, "figures/thesis_gene_expression_enrichment.png"),
       height = 14, width = 30, units = "cm")

5.2.0.1 Previous (discarded) enrichment approaches

# perform_enrichment_analysis <- function(celltype, virustype, target_genes, backgrounds, aquino_data, n_iterations = 1000, seed = 613) {
# 
#   # Select the subset of aquino_data for the given celltype and virustype
#   aquino_subset <- aquino_data %>%
#     dplyr::filter(Cell == celltype, virus == virustype) %>%
#     dplyr::select(gene = GeneID, value = TPM)
# 
#   # Number of target genes
#   n_genes_tested <- length(target_genes)
# 
#   # Identify expressed target genes and number expressed
#   expressed_target_genes <- intersect(target_genes, aquino_subset$gene)
#   n_target_genes_expressed <- length(expressed_target_genes)
# 
#   # Total number of expressed genes in aquino_subset
#   expressed_genes <- aquino_subset$gene
#   M <- length(expressed_genes)
# 
#   # Initialize list to store results for each background
#   results_list <- list()
# 
#   # Loop over each background
#   for (bg_name in names(backgrounds)) {
#    # bg_name <- names(backgrounds)[1]
#     gene_background <- backgrounds[[bg_name]]
# 
#     # Background genes excluding target genes
#     background_only_genes <- setdiff(gene_background, target_genes)
#     N <- length(gene_background)
# 
#     # Set seed for reproducibility
#     set.seed(seed)
# 
#     # Perform random sampling using replicate for efficiency
#     random_counts <- replicate(n_iterations, {
#       # Randomly sample genes from the background
#       random_sample <- sample(background_only_genes, size = n_genes_tested, replace = FALSE)
# 
#       # Count the number of sampled genes that are expressed
#       sum(random_sample %in% expressed_genes)
#     })
# 
#     # Median number of expressed genes in random samples
#     median_random_expressed <- median(random_counts)
# 
#     # Standard deviation
#     sd_random_expressed <- sd(random_counts)
# 
#     # 95% Confidence Interval
#     ci_lower <- quantile(random_counts, 0.025)
#     ci_upper <- quantile(random_counts, 0.975)
# 
#     # Empirical p-value
#     p_value_empirical <- sum(random_counts >= n_target_genes_expressed) / n_iterations
# 
#     # Hypergeometric test p-value
#     n <- n_genes_tested
#     k <- n_target_genes_expressed
#     p_value_hyper <- phyper(k - 1, M, N - M, n, lower.tail = FALSE)
# 
#     # Compile results into a dataframe
#     result_df <- data.frame(
#         label = str_c(celltype, "_", virustype),
#         celltype = celltype,
#         virustype = virustype,
#         background = bg_name,
#         n_genes_tested = n_genes_tested,
#         n_target_genes_expressed = n_target_genes_expressed,
#         median_random_expressed = median_random_expressed,
#         sd_random_expressed = sd_random_expressed,
#         ci_lower = as.numeric(ci_lower),
#         ci_upper = as.numeric(ci_upper),
#         p_value_empirical = p_value_empirical,
#         p_value_hyper = p_value_hyper
#     )
# 
#     results_list[[bg_name]] <- result_df
#   }
# 
#   # Combine the results into one data frame
#   results_df <- do.call(rbind, results_list)
# 
#   return(results_df)
# }

# 
# # Initialize a list to store all results
# all_results <- list()
# 
# # Loop over cell types and virus types
# for (cell in celltypes) {
#   for (virus in virustypes) {
#     # Perform the enrichment analysis
#     result <- perform_enrichment_analysis(
#       celltype = cell,
#       virustype = virus,
#       target_genes = NDD_geneID,
#       backgrounds = enrichment_backgrounds,
#       aquino_data = NDD_data_enrichment,
#       n_iterations = 1000,
#       seed = 123
#     )
# 
#     # Append the result to the list
#     key <- paste(cell, virus, sep = "_")
#     all_results[[key]] <- result
#   }
# }
# 
# # Combine all results into one data frame
# final_results_df <- do.call(rbind, all_results)
# 
# # View the final results
# 
# 
# names(final_results_df)

# There are 12676 unique genes in the Aquino dataset
# long_data$GeneID %>%  unique() %>% length()

# 
# NDD_data_enrichment <- long_data %>% 
#     dplyr::select(GeneID, Cell, virus, TPM)
# 



# Start loop
# 
# 
# celltype = "MONO CD14"
# virustype = "NS"
# gene_background = enrichment_backgrounds[['all_gtex_genes']]
# 
# # Select out df per Cell and virus
# aquino_subset <- NDD_data_enrichment %>% 
#     dplyr::filter(Cell == celltype,
#                   virus == virustype) %>% 
#     dplyr::select(gene = GeneID, value = TPM)
# 
# 
# # Number of target genes
# n_genes_tested <- length(NDD_geneID)
# 
# # Identify expressed target genes and number expressed
# expressed_target_genes <- intersect(NDD_geneID, aquino_subset$gene)
# n_target_genes_expressed <- length(expressed_target_genes)
# 
# # Background genes excluding target genes
# background_only_genes <- setdiff(gene_background, NDD_geneID)
# 
# # Number of iterations for random sampling
# n_iterations <- 1000
# 
# # Initialize vector to store counts of expressed genes in random samples
# random_counts <- numeric(n_iterations)
# 
# # Set seed for reproducibility
# set.seed(613)
# 
# for (i in 1:n_iterations) {
#   # Randomly sample genes from the background
#   random_sample <- sample(background_only_genes, size = n_genes_tested, replace = FALSE)
#   
#   # Count the number of sampled genes that are expressed
#   n_expressed_in_sample <- sum(random_sample %in% aquino_subset$gene)
#   
#   # Store the count
#   random_counts[i] <- n_expressed_in_sample
# }
# 
# 
# # Median number of expressed genes in random samples
# median_random_expressed <- median(random_counts)
# 
# # Standard deviation
# sd_random_expressed <- sd(random_counts)
# 
# # 95% Confidence Interval
# ci_lower <- quantile(random_counts, 0.025)
# ci_upper <- quantile(random_counts, 0.975)
# 
# # Empirical p-value
# p_value_empirical <- sum(random_counts >= n_target_genes_expressed) / n_iterations
# 
# 
# # Hypergeometric Test
# # Total number of background genes
# N <- length(gene_background)
# # Total number of expressed genes in the background
# M <- length(aquino_subset$gene)
# # Number of genes in the sample (same as number of target genes)
# n <- n_genes_tested
# # Number of expressed genes in the target set
# k <- n_target_genes_expressed
# # Hypergeometric test p-value
# p_value_hyper <- phyper(k - 1, M, N - M, n, lower.tail = FALSE)
# 
# # Compile results into a dataframe
# output_df <- data.frame(
#   n_genes_tested = n_genes_tested,
#   n_target_genes_expressed = n_target_genes_expressed,
#   median_random_expressed = median_random_expressed,
#   sd_random_expressed = sd_random_expressed,
#   ci_lower = ci_lower,
#   ci_upper = ci_upper,
#   p_value_empirical = p_value_empirical,
#   p_value_hyper = p_value_hyper
# )
# 
# # Display the results
# print(output_df)
# 
# # Prepare data for plotting
# plot_df <- data.frame(
#   category = c("Target Genes", "Random Samples"),
#   expressed_genes = c(n_target_genes_expressed, median_random_expressed),
#   ci_lower = c(n_target_genes_expressed, ci_lower),
#   ci_upper = c(n_target_genes_expressed, ci_upper)
# )
# 



# 
# # Create dataframe including unexpressed genes (adding rows where TPM=0)
# # Add in background genes and NDD genes that aren't expressed in the enrichment subset (assigning TPM = 0)
# unexpressed_genes <- c(setdiff(gene_background, aquino_subset$gene), # Find genes in background not in aquino subset
#                        setdiff(NDD_geneID, aquino_subset$gene)) %>%  # Find NDD genes not in aquino subset
#     unique()
# 
# 
# ## Background genes
# unexpressed_genes <- tibble(
#     gene = unexpressed_genes,
#     value = 0
# )
# 
# enrichment_subset <- bind_rows(aquino_subset,
#                                unexpressed_genes)
# 
# target_genes <- enrichment_subset %>%
#     dplyr::filter(gene %in% NDD_geneID)
# background_genes <- enrichment_subset %>%
#     dplyr::filter(!gene %in% NDD_geneID)
# 
# 
# # Check length of dataframe
# stopifnot(dim(target_genes)[1] == length(NDD_geneID))
# 
# 
# # Number of expressed genes in the target set
# n_target_genes_expressed <- sum(target_genes$value > 0)
# n_genes_tested <- nrow(target_genes)
# 
# # Number of iterations for random sampling
# n_iterations <- 1000
# set.seed(613)
# 
# # Initialize vector to store counts of expressed genes in random samples
# random_counts <- numeric(n_iterations)
# 
# for (i in 1:n_iterations) {
#   # Randomly sample genes from the background (excluding target genes)
#   random_sample <- background_genes %>%
#     sample_n(size = n_genes_tested, replace = FALSE)
#   
#   # Count the number of expressed genes in the random sample
#   random_counts[i] <- sum(random_sample$value > 0)
# }
# 
# # Median number of expressed genes in random samples
# median_random_expressed <- median(random_counts)
# 
# # Standard deviation of the number of expressed genes in random samples
# sd_random_expressed <- sd(random_counts)
# 
# # 95% Confidence Interval from the random samples (quantiles are empirical confidence intervals in a bootstrapping situation)
# ci_lower <- quantile(random_counts, 0.025)
# ci_upper <- quantile(random_counts, 0.975)
# 
# # Calculate the p-value (proportion of random samples with as many or more expressed genes)
# p_value_empirical <- sum(random_counts >= n_target_genes_expressed) / n_iterations
# 
# # Create the output data frame
# output_df <- data.frame(
#   n_genes_tested = n_genes_tested,
#   n_target_genes_expressed = n_target_genes_expressed,
#   median_random_expressed = median_random_expressed,
#   sd_random_expressed = sd_random_expressed,
#   ci_lower = ci_lower,
#   ci_upper = ci_upper,
#   p_value_empirical = p_value_empirical#,
# #  p_value_hyper = p_value_hyper
# )
# 
# 
# ggplot(plot_df, aes(x = category, y = expressed_genes)) +
#   geom_point(size = 3) +
#   geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
#   ylab("Number of Expressed Genes") +
#   xlab("") +
#   ggtitle("Comparison of Expressed Genes in Target vs. Random Samples") +
#   theme_minimal()

# test_list_for_enrichment = function(target_genes, background_genes, n_iterations=1000, p_cutoff=0.05, set_seed=F, gen_plot=T){
#     #' @param target_genes: Your target genes in `data.frame` form. The function expects two columns, the first `gene` with the gene name, and the second `value` with the value i.e. correlation stat, TPM etc.
#     #' @param background_genes: Your background gene set in `data.frame` form. The function expects two columns, the first `gene` with the gene name, and the second `value` with the value i.e. correlation stat, TPM etc.
#     #' @param n_iterations: The number of iterations/ 'random gene sets'.
#     #' @param p_cutoff: P-value cut-off to call the outcome of the test, default is 0.05.
#     #' @param set_seed: Setting this to TRUE will make your result reproducible i.e it will be the same every time due to identical random gene set selection. Default is FALSE. 
#     #' @param gen_plot: Setting this to TRUE will output a density plot of the random list medians with the target list median indicated on the distribution  
#     #'
#     #' @return Returns results table
#     #'
#     #' @examples
#     #' # implement function with test data - target med > random simulation 
#     #' test_list_for_enrichment(target_genes = data.frame(gene = paste0("ENSG", c(1:100)),
#     #'                                                  value = runif(100, min=0.2, max=1)), 
#     #'                        background_genes = data.frame(gene = paste0("ENSG", c(101:1100)),
#     #'                                                      value = runif(1000, min=0, max=1)), 
#     #'                        n_iterations = 1000, 
#     #'                        p_cutoff = 0.05, 
#     #'                        set_seed = F, 
#     #'                        gen_plot = T)
#     #' 
#     # import libs
#     library(dplyr)
#     library(ggplot2)
#     
#     # if user wants to set seed, then set at 999
#     if(set_seed==T){
#         set.seed(999)
#     }
#     
#     ## check user input ##
#     
#     # ensure target and background list are character vectors
#     stopifnot(is.data.frame(target_genes)==T)
#     stopifnot(is.data.frame(background_genes)==T)
#     
#     # if list is <1, throw error
#     stopifnot(length(target_genes$gene %>% unique()) > 1)
#     
#     # ensure input gene tables contain only unique values
#     target_genes = target_genes %>% 
#         dplyr::distinct()
#     
#     background_genes = background_genes %>% 
#         dplyr::distinct()
#     
#     background_genes_label = "result" 
#     
#     ############################
#     
#     # define rownames for df
#     rn = c('target_set_med', 'random_set_med', 'num_above_rand_med', 'num_below_rand_med', 'p', 'outcome')
#     # generate empty df to hold results
#     res = as.data.frame(matrix(0L, nrow=length(rn), ncol=length(background_genes_label)))
#     # assign row and colnames to empty df
#     colnames(res) = background_genes_label
#     rownames(res) = rn
#     
#     # get median value for target genes
#     pop_med = median(target_genes$value, na.rm=T)
#     
#     # define empty array to hold medians for n=n_iterations random iterations
#     med_arr = replicate(n=as.integer(n_iterations), 0)
#     
#     # define empty matrix to hold full vectors of random gene values for each iteration
#     rand_arr = matrix(0L, nrow=n_iterations, ncol=length(target_genes$gene))
#     
#     # loop over n n_iterations, picking random gene sets n=n_iterations times, and calculating the median value for each
#     # repeat random picking and median computation for n n_iterations
#     for(j in 0:n_iterations){
#         
#         # select random nuclear genes of size of target gene list
#         rand_selection = sample(x=background_genes$gene, size=length(target_genes$gene), replace=F)
#         
#         # get the values for the random genes selected 
#         rand_vals = background_genes %>% 
#             dplyr::filter(gene %in% rand_selection) %>% 
#             dplyr::pull(value)
#         
#         # collect vector of random values for gene set, and add to rand_arr
#         rand_arr[j,] = rand_vals
#         
#         # calculate median of random gene set, and add to med_arr
#         med_arr[j] = median(rand_vals)
#     }
#     
#     # median value of random medians
#     med_of_med_arr = median(med_arr)
#     
#     # generate density plot to show random list mediation distribution
#     if(gen_plot==T){
#         
#         p = data.frame(id=rep("med_arr", length(med_arr)), value=med_arr) %>% 
#             ggplot(data=., aes(x=value, fill=id)) + 
#             theme_bw() +
#             geom_density(alpha=0.75) + 
#             scale_fill_manual(values = my_palette[1]) +
#             theme(legend.position = "none", 
#                   plot.title = element_text(hjust = 0.5, size=11), 
#                   plot.subtitle = element_text(hjust = 0.5, size=8)) + 
#             geom_vline(xintercept=pop_med, linetype="longdash") + 
#             geom_vline(xintercept=med_of_med_arr, linetype="longdash", alpha=0.75, colour='red') +
#             labs(title="Distribution of random gene set medians", subtitle="Dashed black line = target set median\nDashed red line = median of random set medians") + 
#             xlab("Gene set median") + 
#             ylab("Density")
#         
#         plot = print(p)
#     } 
#     
#     # calculate where target list median lies in pop. of random list medians
#     above_pop_med = sum(med_arr > pop_med)
#     below_pop_med = sum(med_arr < pop_med)
#     
#     # assign values to df
#     res['target_set_med', background_genes_label] = pop_med
#     res['random_set_med', background_genes_label] = med_of_med_arr
#     res['num_above_rand_med', background_genes_label] = above_pop_med
#     res['num_below_rand_med', background_genes_label] = below_pop_med
#     
#     # calculate P-value lower tail - if target median is lower than median of random medians, then calculate lower-tail P-value 
#     if(pop_med < med_of_med_arr){
#         p = as.numeric((n_iterations-above_pop_med)/n_iterations)
#         res['p', background_genes_label] = p
#         if(p<=p_cutoff){
#             res['outcome', background_genes_label] = "Target median sig. < than random gene set medians"
#         } else{
#             res['outcome', background_genes_label] = "Target median not sig. diff to random gene set medians"
#         }
#     }
#     
#     # calculate P-value upper tail - if target median is higher than median of random medians, then calculate upper-tail P-value 
#     if(pop_med > med_of_med_arr){
#         p = as.numeric((n_iterations-below_pop_med)/n_iterations)
#         res['p', background_genes_label] = p
#         if(p<=p_cutoff){
#             res['outcome', background_genes_label] = "Target median sig. > than random gene set medians"
#         } else{
#             res['outcome', background_genes_label] = "Target median not sig. diff to random gene set medians"
#         }
#     }
#     
#     res['p_cutoff', background_genes_label] = p_cutoff
#     
#     enrichment_summary <- list("Enrichment results" = res,
#                                "Enrichment plot" = plot)
#     
#     
# }

# enrichment_result <- test_list_for_enrichment(target_genes, 
#                          background_genes,
#                          n_iterations = 1000,
#                          p_cutoff = 0.05,
#                          set_seed = T,
#                          gen_plot = F)
# 

5.3 Explore specificity of gene expression

To identify the specificity of gene expression by cell type and activation state, a specificity metric tau is used. The metric ranges from 0-1, with a higher value indicating greater specificity.

5.3.1 Specificity by cell type

# Tau across treatment
max_proportion_per_treatment_data <- as_tibble( max_proportion_per_treatment_data) %>% 
   mutate(Cell = as.character(Cell))


# Tau by celltypes
wide_data <- max_proportion_per_treatment_data %>%
  dplyr::select(-virus, -log2TPM, -log_mono_nk_ratio, -cd14_mono_prop, -cd56dim_nk_prop, -TPM, -Lineage) %>% 
  pivot_wider(names_from = Cell, values_from = proportion_expression_per_gene)

# Tau by celltypes
wide_data <- max_proportion_per_treatment_data %>%
  dplyr::select(-virus, -log2TPM, -log_mono_nk_ratio, -cd14_mono_prop, -cd56dim_nk_prop, -TPM, -Lineage) %>% 
  tidyr::pivot_wider(names_from = Cell, values_from = proportion_expression_per_gene)

tau_results <- wide_data %>%
  dplyr::rowwise() %>%
  dplyr::mutate(
    tau = {
      x <- dplyr::c_across(where(is.numeric))  # This ensures only numeric columns, typically the cell data, are used.
      x_hat <- x / max(x, na.rm = TRUE)
      tau_value <- sum(1 - x_hat, na.rm = TRUE) / (length(x) - 1)
      tau_value
    }
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::select(GeneID, Symbol, tau_per_cell = tau)
  
# Finding the cell with the maximum proportion_expression_per_gene for each gene
max_cell_per_gene <- max_proportion_per_treatment_data %>%
  as_tibble() %>% 
  dplyr::group_by(GeneID, Symbol) %>%  # Group by GeneID and Symbol to handle multiple genes
  dplyr::summarise(
    Max_Cell = Cell[which.max(proportion_expression_per_gene)]  # Get the cell with the maximum proportion
  ) %>%
  ungroup()  # Ungroup for further operations if necessary



cell_tau <- left_join(tau_results,
                      max_cell_per_gene,
                      by = c("GeneID", "Symbol")) %>% 
  dplyr::mutate(Lineage = "Cell type") %>% 
  dplyr::left_join(.,
             max_proportion_per_treatment_data %>% 
              dplyr::select(Symbol, disease),
            by = "Symbol") %>% 
  dplyr::mutate(Symbol = fct_relevel(Symbol, cell_factor_levels)) %>% 
  group_by(Symbol) %>% 
  slice(1) %>% 
  ungroup()


cell_tau_plot <- cell_tau %>%
  ggplot(aes(x = "Max Cell, τ", y = fct_rev(Symbol), fill = tau_per_cell)) +
  geom_tile() +  # Tiles for heatmap
  geom_text(aes(label = Max_Cell), color = "white", size = 3, vjust = 0.4) +  # Add labels for Max_Cell
  scale_fill_gradient(low = "grey90", high = "darkslategrey", limits = c(0, 1)) +  # Gradient fill
  facet_nested(disease ~ Lineage, scales = "free", space = "free") +
  labs(x = NULL, y = NULL, fill = "τ") +  # Remove x-axis title, adjust legend title
  theme_minimal(base_family = "Helvetica", base_size = 14) +
  theme(
    axis.title.y = element_blank(),  # Remove y-axis title
    axis.text.y = element_blank(),  # Remove y-axis text if unnecessary
    axis.ticks.y = element_blank(),  # Remove y-axis ticks
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
    panel.spacing = unit(0, "lines"),  # Remove space between panels
    strip.background.x = element_rect(fill = NA, colour = "black", linewidth = 1),  # Black box around facet headings
    plot.margin = unit(c(1, 1, 1, 1), "points"),  # Adjust plot margins
    legend.position = "bottom",
    legend.title = element_text(vjust = 0.7),
    legend.key.size = legend_width
  )

cell_tau %>%  gt() %>% 
  opt_interactive()

5.3.1.1 Find number of instances a cell-type contains the maximal gene count

# Find number of cell-types that contain the maximal gene count
cell_tau %>% 
  count(Max_Cell) %>% 
  dplyr::mutate(percentage = n/34*100) %>% 
  gt() %>% 
  opt_interactive()

5.3.1.2 T test comparing specificity (tau) of myeloid cells (Dendritic and monocytes) to NK cells

dc_mono_values <- cell_tau %>% dplyr::filter(Max_Cell %in% c("pDC", "cDC", "MONO CD14", "MONO CD16")) %>% 
                           dplyr::select(tau_per_cell) %>% unlist()

NK_values <-  cell_tau %>%  dplyr::filter(Max_Cell %in% c("NK CD56brt", "NK CD56dim")) %>%
                           dplyr::select(tau_per_cell) %>% unlist()


# Calculate mean and standard deviation for the first group
dc_mono_stats <- cell_tau %>% 
  dplyr::filter(Max_Cell %in% c("pDC", "cDC", "MONO CD14", "MONO CD16")) %>% 
  summarise(
    mean_tau_per_cell = mean(tau_per_cell, na.rm = TRUE),
    sd_tau_per_cell = sd(tau_per_cell, na.rm = TRUE)
  )


# Calculate mean and standard deviation for the second group
NK_stats <- cell_tau %>%
  filter(Max_Cell %in% c("NK CD56brt", "NK CD56dim")) %>%
  summarise(
    mean_tau_per_cell = mean(tau_per_cell, na.rm = TRUE),
    sd_tau_per_cell = sd(tau_per_cell, na.rm = TRUE)
  )

# Compare the two groups using a t-test
t_test_results <- t.test(x = dc_mono_values,
                         y = NK_values,
                         alternative = "two.sided",
                         var.equal = F)
                                                      
t_test_results

## 
##  Welch Two Sample t-test
## 
## data:  dc_mono_values and NK_values
## t = 2.9989, df = 12.77, p-value = 0.01044
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.05330697 0.32980859
## sample estimates:
## mean of x mean of y 
## 0.8234846 0.6319268

5.3.2 Specificity by cell activation state (viral stimulation)

# Tau by activation state

# Find the cell type with the maximum proportion for each gene
max_cell_per_gene <- filtered_data %>%
  group_by(GeneID, Symbol) %>%
  summarise(Max_Cell = Cell[which.max(proportion_expression_per_gene)],
            Max_Proportion = max(proportion_expression_per_gene),
            .groups = 'drop')

# Merge to find corresponding activation state
max_activationState_per_gene <- filtered_data %>%
  inner_join(.,
             max_cell_per_gene, 
             by = c("GeneID", "Symbol", "Cell" = "Max_Cell")) %>%
  dplyr::filter(proportion_expression_per_gene == Max_Proportion) %>%
  dplyr::select(GeneID, Symbol, Cell, virus) %>%
  dplyr::rename(Max_activationState = virus)



# Filter original data to only include entries for the max cell per gene
filtered_max_cells <- filtered_data %>%
  inner_join(max_cell_per_gene, by = c("GeneID", "Symbol", "Cell" = "Max_Cell"))

# Calculate tau for each activationState of these specific max cells
tau_activationStates <- filtered_max_cells %>%
  group_by(GeneID, Symbol) %>%
  summarise(
    tau_activationState = {
      # Normalize the proportions by the maximum observed proportion in the group
      x <- proportion_expression_per_gene
      x_hat <- x / max(x)
      tau_value <- sum(1 - x_hat) / (length(x) - 1)
      tau_value
    },
    .groups = 'drop'
  )


activationState_tau <- left_join(max_activationState_per_gene,
                           tau_activationStates,
                           by = c("GeneID", "Symbol")) %>% 
    dplyr::mutate(Lineage = "Activation state") %>% 
  left_join(.,
            max_proportion_per_treatment_data %>% 
              dplyr::select(Symbol, disease),
            by = "Symbol") %>% 
  dplyr::mutate(Symbol = fct_relevel(Symbol, cell_factor_levels)) %>% 
    group_by(Symbol) %>% 
  slice(1) %>% 
  ungroup() %>% 
  arrange(desc(tau_activationState))

activationState_tau %>% gt() %>% 
  opt_interactive(use_filters = T )

# Create a new column with mapped labels
activationState_tau <- activationState_tau %>%
  mutate(Max_activationState_Label = recode(Max_activationState, !!!plot_labels))

activationState_tau_plot <- activationState_tau %>%
  ggplot(aes(x = "Max Activation state, τ", y = fct_rev(Symbol), fill = tau_activationState)) +
  geom_tile() +  # Tiles for heatmap
  geom_text(aes(label = Max_activationState_Label), color = "white", size = 3, vjust = 0.4) +  # Add labels for Max_activationState
  scale_fill_gradient(low = "grey90", high = "darkslategrey", limits = c(0,1) ) +  # Gradient fill
  facet_nested(disease ~ Lineage, scales = "free", space = "free") +
  labs(x = "Expression by Activation state", y = NULL, fill = "τ") +  # Labels
  theme_minimal(base_family = "Helvetica",
           base_size = 14) +
  theme(
    axis.title = element_blank(),  # Remove y-axis title
      axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
    axis.ticks.y = element_blank(),  # Remove y-axis ticks
    panel.spacing = unit(0, "lines"),  # Remove space between panels
    strip.background.x  = element_rect(fill = NA, colour = "black", linewidth = 1),  # Black box around facet headings
    plot.margin = unit(c(1, 1, 1, 1), "points"), # Adjust plot margins
    legend.position = "bottom",
    legend.title = element_text(vjust = 0.7),
    legend.key.size = legend_width,
     axis.text.y = element_text(hjust = 0.5)  # Center-align y-axis text labels
  )

5.4 Combined plot including gene expression, cell type and activation state.

expression_tau_plot <- (aquino_cell_expression_plot_max  + labs(title = NULL, subtitle = NULL) | 
    cell_tau_plot + theme(strip.text.y = element_blank()) | 
    activationState_tau_plot + theme(axis.text.y = element_text(),
                               )) +
  plot_layout(widths = c(8,1,1), guides = "collect") &
  theme(legend.position = "bottom")

expression_tau_plot

5.5 Differential gene expression of European relative to African populations

table_4b <- read_delim(str_c(args$aquino_eqtl_dir_path, "/aquino_supp_table_S4b.csv"), 
                       skip = 7,
                       col_names = T)  %>% 
 dplyr::filter(Symbol %in% summary_gene_list,
              FDR_lineage_adj <=0.05) %>% 
       #         popDE_lineage_adj == "yes") %>% 
  dplyr::select(-beta_raw, -FDR_raw, -popDE_raw, -t_delta_beta, -popDE_lineage_adj) %>% 
  left_join(.,
            summary_genes,
            by = c("Symbol" = "gene")) %>% 
    relocate(Lineage, disease) 


table_4b_condensed <- table_4b %>% 
  group_by(Lineage, Symbol) %>%
  summarise(
    max_beta_lineage_adj_per_Condition = beta_lineage_adj[which.max(abs(beta_lineage_adj))],
    disease = disease[which.max(abs(beta_lineage_adj))],
    Condition = Condition[which.max(abs(beta_lineage_adj))],
    .groups = 'drop'  # Drops the grouping structure after summarizing
  )


euro_vs_african_expression_plot <- table_4b_condensed %>% 
   ggplot(aes(x = fct_relevel(Lineage, aquino_cell_order$Lineage), 
              y = fct_rev(Symbol), 
              fill = max_beta_lineage_adj_per_Condition)) +
  geom_tile() +
  facet_grid(disease ~ .,
             scales = "free",
             space = "free") +
    scale_fill_gradient2(low = "darkblue", high = "darkred", mid = "white",
                       midpoint = 0,  # This sets the midpoint at 0
                       limit = c(min(table_4b$beta_lineage_adj, na.rm = TRUE), 
                                 max(table_4b$beta_lineage_adj, na.rm = TRUE)),
                       space = "Lab",
                       name = "Log2 fold change") +
    scale_x_discrete(labels = plot_labeller) +
  theme_aw +
  theme(legend.key.size = legend_width) +
   labs(fill = "Log2fold change",
       x = NULL, #"Cell lineage",
       y = "Gene")


euro_vs_african_expression_plot +
  labs(title = "DGE of cell lineages by ancestry", 
       subtitle = "European relative to African")

6 Enrichment of gene expression results

6.1 Enrichment of Chen associations

7 Compiled Figure 1

fig1_gene_expression_chen <- 
  (
    (aquino_gene_count + theme_aw + theme(legend.direction = "vertical", 
                                          legend.title = element_text(vjust = 1)) | 
       chen_gene_count +theme_aw + theme(legend.position = "bottom",
                                         legend.title = element_text(vjust = 1)) | 
       euro_vs_african_expression_plot + theme_aw + theme(legend.position = "bottom",
                                                          panel.grid.major = element_blank(),
                                                          panel.grid.minor = element_blank(),
                                                          legend.key.width = unit(0.8, "cm"),
                                                           legend.title = element_text(vjust = 1)) 
    )+ 
      plot_layout(guides = "keep", widths = c(1,4,1.5))
  ) /
  (
    (aquino_cell_expression_plot_max  + theme_aw + 
       theme(legend.key.width = unit(0.8, "cm"),
             legend.title = element_text(vjust = 1))  + 
       labs(title = NULL, subtitle = NULL) | 
       cell_tau_plot + theme_aw +theme(strip.text.y = element_blank(),
                                       legend.key.width = unit(0.8, "cm"),
                                       legend.title = element_text(vjust = 1)) | 
       activationState_tau_plot + labs(x = NULL) + theme_aw +  theme(axis.text.y = element_blank(),
                                                               axis.ticks.y = element_blank(), 
                                                               legend.key.width = unit(0.8, "cm"),
                                                               legend.title = element_text(vjust = 1)
                                                               )
    ) +
      plot_layout(widths = c(10,1.6,1.6), guides = "collect") 
  ) +
  plot_layout(heights = c(3,8)) +
  plot_annotation(tag_levels = list(c("a", "b", "c", "d", "e", "f")
  )
  ) 


fig1_gene_expression_chen

ggsave(fig1_gene_expression_chen,
       height = 36, width = 32, units = "cm",
       file = file.path(base_dir, "figures/fig_1_gene_expression_chen.pdf"))

ggsave(fig1_gene_expression_chen,
       height = 36, width = 32, units = "cm",
       file = file.path(base_dir, "figures/fig_1_gene_expression_chen.png"))

8 Individual figure panels

thesis_fig_3 <- (aquino_gene_count + theme_aw + theme(legend.direction = "vertical",
                                                      legend.title = element_text(vjust = 1)) | 
                     enrichment_plot) +
    plot_annotation(tag_levels = "a")

ggsave(thesis_fig_3,
       file = file.path(base_dir, "figures/thesis3_aquino_gene_count.pdf"),
       width = 18, height = 11, units = "cm")

thesis_figure_chen <- (chen_gene_count + 
                           theme_aw + 
                           theme(legend.direction = "vertical",          
                                 legend.position = "right",
                                 legend.title = element_text(vjust = 1))) /
         chen_association_count +
    plot_annotation(tag_levels = "a") +
    plot_layout(heights = c(1,4))

ggsave(thesis_figure_chen,
       file = file.path(base_dir, "figures/thesis_chen_associations.png"),
       width = 38, height = 45, units = "cm")

thesis_figure_chen

thesis_expression_tau_plot <-  (aquino_cell_expression_plot_max  + theme_aw + 
       theme(legend.key.width = unit(0.8, "cm"),
             legend.title = element_text(vjust = 1))  + 
       labs(title = NULL, subtitle = NULL) | 
       cell_tau_plot + theme_aw +theme(strip.text.y = element_blank(),
                                       legend.key.width = unit(0.8, "cm"),
                                       legend.title = element_text(vjust = 1)) | 
       activationState_tau_plot + labs(x = NULL) + theme_aw +  theme(axis.text.y = element_blank(),
                                                               axis.ticks.y = element_blank(), 
                                                               legend.key.width = unit(0.8, "cm"),
                                                               legend.title = element_text(vjust = 1))
                                                               
 ) + plot_layout(widths = c(8,1,1), guides = "collect") +
     plot_annotation(tag_levels = "a") &
  theme(legend.position = "bottom")

thesis_expression_tau_plot

ggsave(thesis_expression_tau_plot,
       file  = file.path(base_dir, "figures/thesis_expression_tau.png"),
       units = "cm", width = 36, height = 30)

ggsave(euro_vs_african_expression_plot,
       file = file.path(base_dir, "figures/thesis_eur_afr_dge.png"),
       units = "cm", width = 12, height = 16)

sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 14.6.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] httr_1.4.7        gt_0.11.0         forcats_1.0.0     stringr_1.5.1    
##  [5] dplyr_1.1.0       purrr_1.0.1       readr_2.1.5       tidyr_1.3.1      
##  [9] tibble_3.2.0      tidyverse_1.3.2   scales_1.3.0      readxl_1.4.1     
## [13] ggh4x_0.2.8       ggplot2_3.4.2     viridis_0.6.2     viridisLite_0.4.1
## [17] wesanderson_0.3.6 paletteer_1.4.1   sciRmdTheme_0.1   patchwork_1.2.0  
## [21] here_1.0.1       
## 
## loaded via a namespace (and not attached):
##  [1] fs_1.6.4            reactable_0.4.4     lubridate_1.8.0    
##  [4] bit64_4.0.5         rprojroot_2.0.3     tools_4.2.0        
##  [7] backports_1.4.1     bslib_0.3.1         utf8_1.2.2         
## [10] R6_2.5.1            DBI_1.2.2           colorspace_2.1-0   
## [13] withr_2.5.0         tidyselect_1.2.1    gridExtra_2.3      
## [16] processx_3.8.4      bit_4.0.5           compiler_4.2.0     
## [19] textshaping_0.3.6   cli_3.6.3           rvest_1.0.4        
## [22] xml2_1.3.6          labeling_0.4.3      bookdown_0.29      
## [25] sass_0.4.1          systemfonts_1.0.6   digest_0.6.35      
## [28] rmarkdown_2.14      pkgconfig_2.0.3     htmltools_0.5.7    
## [31] highr_0.9           dbplyr_2.2.1        fastmap_1.1.1      
## [34] htmlwidgets_1.6.4   rlang_1.1.1         rstudioapi_0.13    
## [37] farver_2.1.1        jquerylib_0.1.4     generics_0.1.3     
## [40] jsonlite_1.8.8      crosstalk_1.2.1     vroom_1.6.5        
## [43] googlesheets4_1.1.1 magrittr_2.0.3      Rcpp_1.0.12        
## [46] munsell_0.5.1       fansi_1.0.3         lifecycle_1.0.3    
## [49] stringi_1.7.6       yaml_2.3.5          grid_4.2.0         
## [52] parallel_4.2.0      promises_1.3.0      crayon_1.5.2       
## [55] haven_2.5.4         chromote_0.3.1      hms_1.1.2          
## [58] knitr_1.40          ps_1.7.6            pillar_1.9.0       
## [61] reprex_2.1.0        glue_1.6.2          evaluate_0.23      
## [64] modelr_0.1.11       vctrs_0.6.5         tzdb_0.4.0         
## [67] cellranger_1.1.0    gtable_0.3.1        reactR_0.6.0       
## [70] rematch2_2.1.2      assertthat_0.2.1    xfun_0.34          
## [73] broom_1.0.5         later_1.3.2         ragg_1.2.5         
## [76] googledrive_2.1.0   gargle_1.5.0        websocket_1.4.1    
## [79] ellipsis_0.3.2